Earth Stress Management and Control Process For Hydrocarbon Recovery

ABSTRACT

A method for controlling hydrocarbon recovery to improve well interactions while preventing excessive strain or stress-induced well deformations and mechanical failures is described. The method incorporates a systematic, transient analysis process for determining the formation effective displacement, stress and excess pore pressure field quantities at any depth within a stratified subterranean formation resulting from the subsurface injection and/or withdrawal of pressurized fluids.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Application No.60/845,951, filed 20 Sep. 2006.

U.S. provisional patent application, attorney docket number 2006EM115,entitled “Fluid Injection Management Method for Hydrocarbon Recovery,”and U.S. provisional patent application, attorney docket number2006EM117, entitled “Earth Stress Analysis Method for HydrocarbonRecovery,” filed concurrently herewith, contain subject matter relatedto that disclosed herein, and are incorporated herein by reference intheir entirety.

BACKGROUND OF THE INVENTION

1. Field of the Invention

Embodiments of the present invention relate generally to the analysis ofearth stresses associated with hydrocarbon recovery and, moreparticularly, to determining effective displacements and stressesresulting from the injection and/or withdrawal of pressurized fluids.

2. Description of the Related Art

Hydrocarbon recovery processes may occur in subterranean reservoir sandsor shales, may employ single or multiple water injection wells, and maybe energetic by design. Recovery processes may require large pressureand/or temperature changes to promote the extraction of oil and gas ateconomic rates from subterranean formations. Water injection wells maybe employed in either a secondary water flood strategy to sweep residualoil to production wells and provide pressure support, or strictly todispose of produced water. The pressure or temperature changes, inducedduring the injection and/or withdrawal of pressurized fluids, result instresses to the formation that may lead to some combination offracturing, expansion and contraction of the subterranean formations.

These fracturing, expansion and contraction processes typically causeexcess pore pressure and stress changes within the formations that maybe large enough to negatively impact well mechanical integrity,productivity, injectivity and conformance. They may also be large enoughto exceed the mechanical limits of penetrating wells. If the mechanicallimits are exceeded, some combination of expansion and fracturing of thewell or subterranean formation may occur. As a result, the penetratingwells may no longer be capable of sustaining reliable hydrocarbonproduction safely and without risk to the environment.

Many of the same risks are present in water flood or water disposalcampaigns and success may depend on the capabilities to manage earlywater breakthrough and contain hydrofracture growth within the targetsubterranean interval(s). For water flooding, the expansion andfracturing process may lead to “short circuiting” of injector-producerpatterns and loss of pressure support. For water disposal, theseprocesses may lead to a loss of containment that may result inrepressurization of untargeted zones and, potentially, regulatory andenvironmental consequences.

Prior art methods employed for analysis of earth stresses associatedwith the injection or production of hydrocarbons have usually adoptedone of two approaches. In the first approach, conventional well logging(e.g., gamma ray, density, resistivity, and sonic) analysis techniquesare utilized in conjunction with production data to infer changes inearth stresses. In the second approach, earth stresses are determined byanalytic models or simulators. Either approach typically assumes steadystate conditions and is specific only to a particular set of wellperformance conditions (e.g., single-valued average pressure, rate,temperature and single-layered formation properties).

These conventional approaches fall short of being generalized to accountfor multiple subterranean layers and variable, time-dependent wellperformance. In addition, even though displacement measurements fromfield surveillance may indicate the presence of multi-well interactions,conventional approaches do not scale very easily to account for theseinteractions at the field-level. Field surveillance methods may includesurveys of ground surface displacements via tilt arrays, remote sensing(e.g., Interferometric Synthetic Aperture Radar (InSAR), Light Detectionand Ranging (LiDAR), Global Positioning System (GPS)) or verticalprofiling of recorded passive and/or active microseismic (μ-seismic)events. The underlying methodology of the prior art also precludes rapidforward or inverse modeling with field surveillance data to furtherconstrain modeling problems and allow calibration of the model withcollected field surveillance data.

The prior art detailing the methods of controlling subterraneaninjection and hydrocarbon production processes has not focused onmulti-well control or the enablement of field-wide control systems.Moreover, the scope of the prior art has been limited to detecting sometime-dependent, single-well characteristic of the resident or injectedfluids, changes in the geometry of a hydrofracture, or a principalstress change within the very near-well regime for predicting phenomena,such as the potential for or onset of sand production.

Accordingly, what is needed is a well-based and/or a field-based,injection control process that accurately models multi-layeredsubterranean formations and predicts injection conditions required toimprove injector performance while minimizing undesirable fracturegrowth and the potential for loss of fracture containment.

SUMMARY OF THE INVENTION

One embodiment of the present invention is a method for managing animpact, on an earth formation, of water injection operations associatedwith hydrocarbon recovery from at least one well formed in the earthformation. The method generally includes dividing the well into aplurality of layers; generating at least first and second sets ofequations to model contributions to the impact of the operations due toat least first and second physical processes associated with theoperations; obtaining solutions to the first and second sets ofequations to determine contributions to the impact of the operations dueto the first and second physical processes; combining the solutions tothe first and second sets of equations to determine the impact of theoperations on the earth formation; forecasting an injection mode of thewater injection operations; and using earth surface displacementmeasurements and data gathered while monitoring the water injectionoperations to update the forecasted mode.

Another embodiment of the present invention is a method of determiningan impact, on an earth formation, of water injection operationsassociated with hydrocarbon recovery from a field formed in the earthformation, wherein the field comprises a plurality of wells. The methodgenerally includes dividing each of the wells into a plurality oflayers; generating at least first and second sets of equations for eachof the layers to model contributions to the impact of the operations dueto at least first and second physical processes associated with theoperations; obtaining solutions to the first and second sets ofequations to determine contributions to the impact of the operations dueto the first and second physical processes for each of the layers;combining the solutions to the first and second sets of equations toobtain combined layer solutions for each of the layers; superposing thecombined layer solutions for each of the layers to obtain combined wellsolutions for each of the wells; superposing the combined well solutionsto determine the impact of the operations on the earth formation;forecasting the injection mode of the water injection operations; andusing earth surface displacement measurements and data gathered whilemonitoring the operations to update the forecasted mode.

Yet another embodiment of the present invention is a method for managingan impact, on an earth formation, of fluid injection operationsassociated with hydrocarbon recovery from at least one well formed inthe earth formation. The method generally includes dividing the wellinto a plurality of layers; generating at least first and second sets ofequations to model contributions to the impact of the operations due toat least first and second physical processes associated with theoperations; obtaining solutions to the first and second sets ofequations to determine contributions to the impact of the operations dueto the first and second physical processes; combining the solutions tothe first and second sets of equations to determine the impact of theoperations on the earth formation; forecasting the injection mode of thefluid injection operations; and using earth surface displacementmeasurements and data gathered while monitoring the fluid injectionoperations to update the forecasted mode.

BRIEF DESCRIPTION OF THE DRAWINGS

So that the manner in which the above recited features of the presentinvention can be understood in detail, a more particular description ofthe invention, briefly summarized above, may be had by reference toembodiments, some of which are illustrated in the appended drawings. Itis to be noted, however, that the appended drawings illustrate onlytypical embodiments of this invention and are therefore not to beconsidered limiting of its scope, for the invention may admit to otherequally effective embodiments.

FIG. 1 illustrates a decomposition of a field-scale, multi-well probleminto single-well effects in accordance with an embodiment of the presentinvention.

FIG. 2 is a flowchart illustrating an analysis method for determiningthe formation effective displacement, stress and excess pore pressurefield quantities resulting from injection and/or production inaccordance with an embodiment of the present invention.

FIG. 3 is a flowchart illustrating a systematic method for passivemicroseismic (μ-seismic) profiling in accordance with an embodiment ofthe present invention.

FIG. 4 is a flowchart illustrating integration of the methods in FIGS. 2and 3 in accordance with an embodiment of the present invention.

FIG. 5 depicts a superposed, single-well solution for the displacementand stress field quantities in accordance with an embodiment of thepresent invention.

FIG. 6 illustrates a model predictive control (MPC) scheme to forecastwell operating conditions relative to well mechanical integrity inaccordance with an embodiment of the present invention.

FIG. 7 illustrates an MPC scheme that incorporates earth displacementmeasurements in accordance with an embodiment of the present invention.

FIG. 8 illustrates an MPC scheme that simultaneously incorporates earthdisplacement measurements and vertical profiling of μ-seismic events inaccordance with an embodiment of the present invention.

FIG. 9 illustrates a flow chart of a computer-based MPC scheme for inaccordance with an embodiment of the present invention.

FIG. 10 illustrates a main graphical user interface (GUI) forcharacterization and control of earth stresses associated with ahydrocarbon recovery process in accordance with an embodiment of thepresent invention.

FIG. 11 illustrates a data management module in the main GUI of FIG. 10in accordance with an embodiment of the present invention.

FIG. 12 illustrates an execution form of a data management module in themain GUI of FIG. 10 in accordance with an embodiment of the presentinvention.

FIG. 13 illustrates a well log calculation module in the main GUI ofFIG. 10 in accordance with an embodiment of the present invention.

FIG. 14 illustrates a field quantity calculation module in the main GUIof FIG. 10 in accordance with an embodiment of the present invention.

FIG. 15 illustrates an execution form of a field quantity calculationmodule in the main GUI of FIG. 10 in accordance with an embodiment ofthe present invention.

FIG. 16 illustrates a field view module in the main GUI of FIG. 10 inaccordance with an embodiment of the present invention.

FIG. 17 illustrates a well integrity module in the main GUI of FIG. 10in accordance with an embodiment of the present invention.

FIG. 18 illustrates a shear-slip limit approach used in a well integritymodule in the main GUI of FIG. 10 in accordance with an embodiment ofthe present invention.

FIG. 19 illustrates an MPG self-calibration module in the main GUI ofFIG. 10 in accordance with an embodiment of the present invention.

FIG. 20 illustrates an MPG scheme to forecast well operating conditionsrelative to injection constraints in accordance with an embodiment ofthe present invention.

FIG. 21 illustrates an MPG scheme that incorporates earth displacementmeasurements to self-calibrate in accordance with an embodiment of thepresent invention.

FIG. 22 illustrates an MPG scheme that simultaneously incorporates earthdisplacement measurements and vertical profiling of μ-seismic events inaccordance with an embodiment of the present invention.

DETAILED DESCRIPTION

Embodiments of the present invention provide a systematic, transientanalysis method for determining the formation effective displacement,stress and excess pore pressure field quantities at any depth within astratified subterranean formation resulting from the subsurfaceinjection and/or withdrawal of pressurized fluids; a process forcontrolling recovery to improve well interactions while preventingexcessive strain or stress-induced well deformations and mechanicalfailures; and a process for controlling fluid injection parameters toimprove well interactions and control hydrofracture geometries.

Embodiments of the present invention also incorporate data from fieldsurveillance, well logs, well trajectories, completions and/or variousother injection or production sources for controlling various wellparameters, and in self-calibrating the model. Water flood and/or waterdisposal operations may also be considered for some embodiments increating or evaluating a fieldwide development strategy.

An Exemplary Earth Stress Analysis Method

The disclosed analysis method provides a logical sequence fordetermining said field quantities given a representative set of welllogs, well trajectories, completion types and injection or productiondata (e.g., pressures, rates, temperatures and fluid properties). Indetermining said field quantities, multi-well solution methods may bederived by superposing physics-based single-well solution methods ofgoverning processes, such as poroelastic expansion or contraction,thermoelastic expansion or contraction, and dislocations or fractures.Superposing, in such a case, may be considered as a summation of thecalculations of various forces, initially considered independently, at aspecified location.

The physics-based solution method for analysis of the multi-well problemmay be based on mathematical decomposition of the aforementionedgoverning processes on a single well basis. As an example, the case of asubsurface injection and/or production process which employs a heatedfluid as an energetic injectant to aid in recovering hydrocarbons fromsubterranean formations may be considered. The complex recovery processmay be systematically decomposed into constitutive effects in which thephysics governing the effects may be well understood. As shown in FIG.1, a full field, multi-well interactions model 100 may be decomposedinto multiple single-well representations, and the single-wellrepresentations 110 may be further decomposed into some combination ofvarious effects, such as poroelastic, thermoelastic and dislocationeffects 170, 180, 190.

For the example shown in FIG. 1, the conditions may be predicated uponthe assumptions of energetic injectant pressure 130 and temperatureabove initial subterranean reservoir conditions, axisymmetry about asingle well 120 and in-situ stress conditions favoring the initiation ofhydraulically-induced horizontal fractures 150, but these assumptionsmay not be required for general consideration in determining fieldquantities. In fact, the systematic earth stress analysis method may beflexible enough to account for arbitrary boundary conditions 140,in-situ stress states, injectant conditions and orientations ofdislocations or fractures (if initiated). The method may also beflexible enough to account for the poroelastic 170, thermoelastic 180and dislocation effects 190 singularly, or as a combination of saideffects.

For injection, the induced subterranean formation dilation andfracturing 160 may be decomposed into the equivalent effects ofporoelastic expansion 170, thermoelastic expansion 180 (if the processis thermal), and opening and/or shear dislocation 190 (if fracturingoccurs). For production, fracturing 160 may also be decomposed vis-à-visinjection. In either case, a systematic sequence of calculations may bemade on a single-well basis for the effective components of displacementand stress as shown in FIG. 2. Data may be collected from adepth-dependent log, which may have only sparse information and onlyinclude gamma ray measurements, or the log information may includedensity, resistivity and sonic (monopole or dipole) measurements.

The principal assumptions that may be required for decomposition of theinjection and production problems into constitutive poroelastic andthermoelastic effects may be the following: the injected and producedfluids may be incompressible, may be Newtonian and flow from pointsources. With these assumptions, the injection rates may be treated aspiecewise constant within the smallest time interval, but may bevariable over longer intervals of time. The subterranean earth model maybe composed of multiple, transversely isotropic layers and may be viewedmathematically as a propagation of layered elastic half-space solutions.The layered earth model may be prestressed with a uniform lateralcompressive stress (σ₀) and an axial stress equivalent to weight (ρgh)of the overlying strata.

Fractures may initiate instantaneously when injection pressure risesabove a local fracture gradient and may close instantaneously wheninjection pressure falls below the local fracture gradient. Fractureleakoff and thermal conduction may be normal to local fracture faces,and fracture loading is symmetric without tip effects. The radial extentof pressure and thermal fronts when injecting below a local fracturegradient may be dictated by ordinary diffusion processes. The radialextent of a pressure front when injecting above a local fracturegradient may be considered equivalent to the fracture radius. The radialextent of a thermal front when injecting a heated fluid may beequivalent to the limit of advance within fractures.

Mode I opening of fractures, where the walls of the opening moveperpendicularly away from the fracture plane when the fracture formed,may be equivalent to the normal displacement discontinuity. Mode IIopenings, or openings due to in-plane shear, may be equivalent to theshear displacement discontinuity. Shear stress at the free surface ofthe layered earth model may be zero. The injection-induced fractureproblem may be equivalent to superposition of the poroelastic problem,the thermoelastic problem, and the dislocation problem (see FIG. 1).

FIG. 2 illustrates the general method 200 for analysis of earthstresses. The fundamental analysis may begin with the single-layeredelastic half-space solution given representative input data, andsystematic workflow may be established to decompose the injection orproduction problem into constitutive effects. The illustrated methodbegins at block 202, and inputs may be read at block 204. The inputs mayinclude formation properties, such as rock data 208 (e.g., layer elasticproperties, layer fracture toughness, and the like) and fluid data 206(e.g., energetic fluid properties, fluid rate and pressure). The inputsmay also include cycle data 210 and stimulation data 212. In variousoperations, including those involving employing a pressurized injectant,such as steam, injectant pressure 216 and injectant rate 218 may beread. The method may be designed to analyze the data collected over aspecified period of time for a particular well and then proceed toanalyze the data for another well, looping through the wells in turn. Atblock 220 it may be determined if there is a well to analyze. If thereis no well to analyze or if all wells have already been analyzed, themethod may stop or may proceed to block 222.

If there are more wells to analyze, the method may proceed to block 222it may be determined if there is a time increment of data to analyze. Ifall the time increments of data have been analyzed, then the method mayreturn to block 220 where the next well to be analyzed may be selected.If there is a time increment to analyze, it may be determined whetherthere is a flow rate at block 223. The flow rate may be the flow rate ofsteam. If there is a flow rate, various data, such as an oil flow rate228 and a water flow rate 230, may be read into the system at block 226.At block 224 the pressure may be analyzed. Blocks 232 and 234 in thesequence may be performed in an effort to determine the fracture extent,fracture width and thermal extent at time t if injection is above thelocal fracture gradient. If fracturing and thermal effects are ignoredthen the pressure extent is evaluated using ordinary diffusionrelationships.

For fracturing, the extent may be calculated via a convolution of arepresentative solution (Carter, R. D., Derivation of the GeneralEquation for Estimating the Extent of the Fractured Area, Drill. & Prod.Prac., API (1957); Geertsma, J. & L. R. Kern: Widths of HydraulicFractures, J. Pet. Tech. (September 1961), Trans., AIME) for variablerate, and the convolved width may be calculated in terms of the extentaccordingly. If, for example, the Carter solution is adopted as thesolution for fracture extent, then the corresponding thermal extent maybe determined via a convolution of the Marx-Langenheim solution (Marx,J. W. & R. H. Langenheim: Reservoir Heating by Hot Fluid Injection”,Trans., AIME, Vol. 216 (1959)), which may be made analogous to theCarter solution. At block 236, the pressure and temperature gradientsmay be evaluated starting from time-dependent (preferentially real-time)pressure and temperature measurements.

If both of these measurements are not available for injection, asuitable starting point may then become the isobaric or isothermalsaturation values. For production, the pressure and temperaturegradients may also be evaluated starting from time-dependent(preferentially real-time) pressure and temperature measurements. Ifthese measurements are not available for production, then the startingpoint may be derived from a suitable convolution for the gradients interms of isothermal bulk compressibility and reservoir heat loss due toproduction.

The elastic, half-space solution (single-layered or multi-layered) forthe displacement and stress field quantities may be determined as givenat block 238 and may be evaluated in terms of Lipschitz-Hankel typeintegrals I(a,b;d) or the modified type Ĵ_(mnp) involving Besselfunctions J_(a,b,m,n) given by Equation 1 as follows:

$\begin{matrix}{{{I\left( {a,{b;d}} \right)} = {\int_{0}^{\infty}{^{{- q}\; \alpha}\alpha^{d}{J_{a}\left( {\alpha \; R} \right)}{J_{b}\left( {\alpha \; r} \right)}{\alpha}}}}{{\hat{J}}_{mnp} = {\left\{ {{sign}\left( {ϛ - ϛ^{\prime}} \right)} \right\}^{({m + n + p})}{\int_{0}^{\infty}{{J_{m}\left( {t\; \rho} \right)}{J_{n}(t)}^{{- t}{{ϛ - ϛ^{\prime}}}}t^{p}{t}}}}}} & (1)\end{matrix}$

where r is the radial coordinate and R is the fracture or thermalextent; q is (z±h) and (ζ−ζ′) is given by (z±h)/R, where z is thevertical coordinate and h is the burial depth, and the indices a,b;d orm,n,p are 0 to 2.

Blocks 240-244 in the sequence may be performed in parallel to thepreviously described blocks 232-238 (e.g., a production cycle follows aninjection cycle). Otherwise only blocks 232-238 (for injection) orblocks 240-244 (for production) may be required. In either case, thesolutions for injection and production may still be evaluated in termsof the Lipschitz-Hankel type integrals given by Equation 1. A rigorousmathematical formulation for displacements and stresses at any depthwithin the same stratified subterranean formation due to propagation ofinjection induced fractures below the surface may also be developed.

If a solution, such as that provided for in the previously describedCarter reference, is adopted to determine rate and time-dependentfracture extent based on the input data, then a convolution is requiredsince the solution is formulated on the basis of constant rate. Thepreferred convolution for the fracture extent is then given by Equation2 as follows:

$\begin{matrix}{{{\pi \; {R_{F}^{2}\left( t_{n} \right)}} = {{A_{F}\left( t_{n - 1} \right)} + \begin{bmatrix}{{\left( \frac{Q\left( t_{n} \right)}{\pi \; \kappa \; \Delta \; {P\left( t_{n} \right)}} \right)\sqrt{\left( {\pi \; {Dt}_{n}} \right)}} -} \\{\left( \frac{Q\left( t_{n} \right)}{\pi \; \kappa \; \Delta \; {P\left( t_{n} \right)}} \right)\sqrt{\left( {\pi \; {Dt}_{n - 1}} \right)}}\end{bmatrix}}}{{\kappa = \frac{K}{\mu}},{D = {\frac{K_{v}}{K_{h}}\frac{K}{{\mu\varphi}\; C_{t}}}},{C_{t} = \frac{3\left( {1 - {2v}} \right)}{E}}}} & (2)\end{matrix}$

where R_(F) is the half-length of injection induced thermal extent inmeters, t is time in days, n is the index for time, A_(F) is the area ofinjection induced fracture in square meters, Q is the injection orproduct ion rate in cubic meters per day, ΔP is (P_(inj)−P_(res)) or(P_(prod)−P_(res)) in pascals, P_(inj) is the injection pressure inpascals, P_(res) is the initial reservoir pressure in pascals, P_(prod)is production pressure in pascals, κ is the effective mobility to waterin square meters per pascal seconds, D is the pore fluid pressurediffusivity, K is formation permeability in millidarcy, μ is theviscosity of injectant in centipoise, K_(ν)/K_(h) is the permeabilityratio, C_(t), is the formation bulk compressibility, φ is formationporosity in porosity units, E is the formation elastic modulus inpascals, and ν is the formation Poisson's ratio.

Since the fracture width is assumed to be a function of the extent, thewidth solution is naturally rate and time-dependent. For example, thesolution due to the width is given by:

$\begin{matrix}{{{W_{F}\left( {r,t} \right)} = {\frac{\begin{matrix}{8\left( {1 - v^{2}} \right)} \\\sqrt{\left( {R_{F}^{2} - r^{2}} \right)}\end{matrix}}{\pi \; E}\left( {{P_{inj}(t)} - \sigma_{3}} \right)\left( {1 - \frac{A_{pe}}{2}} \right)}},{A_{pe} = \frac{\alpha_{b}\left( {1 - {2v}} \right)}{\left( {1 - v} \right)}}} & (3)\end{matrix}$

where α_(b) is the Biot coefficient.

The solutions to Equations 2 and 3 may be constrained according towhether there is enough pressure available within any time interval toovercome the minimum principal stress local to the point where afracture may initiate and propagate. Therefore, it may be expected thata fracture will initiate and propagate when the criterion given byEquation 4 is satisfied such that

$\begin{matrix}{{{P_{fp} = {P_{foc} + \left\lbrack \frac{\pi \; E\; \gamma_{se}}{2{R_{f}\left( {1 - v^{2}} \right)}} \right\rbrack^{1/2}}},{P_{foc} = {S_{grad}*H}}}{\gamma_{se} = {{\frac{K_{IC}^{2}{\pi \left( {1 - v^{2}} \right)}}{2E}\therefore P_{fp}} = {P_{foc} + {\left\lbrack \frac{\pi^{2}K_{IC}^{2}}{4R_{f}} \right\rbrack^{1/2}.}}}}} & (4)\end{matrix}$

where P_(fp) is the fracture propagation in pressure in pascals, P_(foc)is the opening/closing pressure in pascals, S_(grad) is the maximumprincipal stress gradient in pascals per meter, H is the source burialdepth in meters and K_(IC) is the formation fracture toughness.

Analogous to the solution for rate and time-dependent fracture extent,the rate and time-dependent thermal extent may also be determined. Asolution, such as that provided for in the previously describedMarx-Langenheim reference, is adopted in this example and given byEquation 5, where the temperature profile ahead of the thermal front,but within the fracture, is assumed to be governed by the ordinaryrelationship for thermal conduction in a semi-infinite medium (i.e.,Equation 6).

$\begin{matrix}{{\pi \; R_{T}^{2}} = {{A_{T}\left( t_{n - 1} \right)} + \left\lbrack {{\left( \frac{{Q\left( t_{n} \right)}\rho \; h_{i}}{\pi \; k\; \Delta \; {T\left( t_{n} \right)}} \right)\sqrt{\left( {\pi \; \alpha \; t_{n}} \right)}} - {\left( \frac{{Q\left( t_{n} \right)}\rho \; h_{i}}{\pi \; k\; \Delta \; {T\left( t_{n} \right)}} \right)\sqrt{\left( {\pi \; \alpha \; t_{n - 1}} \right)}}} \right\rbrack}} & (5) \\{{T_{F}\left( {{r \geq R_{T}},t} \right)} = {{\left( {T_{inj} - T_{res}} \right){{erfc}\left( \frac{r}{\sqrt{4\alpha}t} \right)}} + T_{res}}} & (6)\end{matrix}$

where R_(T) is the half-length of injection induced thermal extent inmeters, A_(T) is the area of thermal advancement in square meters, ρ isthe density of injectant in kilograms per meter, h_(i) is the enthalpyof injectant in kilojoules per kilogram, k is the thermal conductivityin watts per meter per degrees Celsius, α is the thermal diffusivity insquare meters per second, and ΔT is (T_(inj)−T_(res)) or(T_(prod)−T_(res)), where, T_(inj), is the injection temperature,T_(res) is the initial reservoir temperature and T_(prod) is theproduction temperature, all in degrees Celsius.

Similarly the gradient and vertical extents of the pressure and thermalfronts, relative to fracture surfaces, may also be evaluated on thebasis of a semi-infinite medium according to the following (Equations 7and 8):

$\begin{matrix}{{{P\left( {z,t} \right)} = {{\left( {P_{inj} - P_{res}} \right){{erfc}\left\lbrack \frac{z}{\left\lbrack {4\; D_{P}t} \right\rbrack^{1/2}} \right\rbrack}} + P_{res}}}{{z_{P}\left( {P,t} \right)} = {2*\left( {{{erfc}^{- 1}\left\lbrack \frac{\left( {P_{res} + P_{inc}} \right) - P_{res}}{\left( {P_{inj} - P_{res}} \right)} \right\rbrack}\left\lbrack {4D_{P}t} \right\rbrack}^{1/2} \right)}}} & (7) \\{{{T\left( {z,t} \right)} = {{\left( {T_{inj} - T_{res}} \right){{erfc}\left\lbrack \frac{z}{\left\lbrack {4D_{T}t} \right\rbrack^{1/2}} \right\rbrack}} + T_{res}}}{{z_{T}\left( {T,t} \right)} = {2*\left( {{{erfc}^{- 1}\left\lbrack \frac{\left( {T_{res} + T_{inc}} \right) - T_{res}}{\left( {T_{inj} - T_{res}} \right)} \right\rbrack}\left\lbrack {4D_{T}t} \right\rbrack}^{1/2} \right)}}} & (8)\end{matrix}$

where D_(p) is the pore fluid pressure diffusivity in square meters persecond and D_(T) is the thermal diffusivity in square meters per second.

If time-dependent temperature data for the injected and/or producedfluid conditions, sampled at reasonably repetitive intervals, is notavailable then it may be plausible to approximate the conditions. Forexample, if steam is the injectant and a constant steam quality isassumed, then the pressure changes associated with injection orproduction may lead to changes in temperature given by the followingEquation 9:

$\begin{matrix}{\begin{matrix}{\left. {\Delta \; {T\left( t_{n} \right)}}\rightarrow{f\left( {\Delta \; {P_{inj}\left( t_{n} \right)}} \right)} \right. = \frac{T_{s}}{T^{\star}}} \\{= \left\lbrack \frac{n_{10} + D - \sqrt{\left( {n_{10} + D} \right)^{2} - {4\left( {n_{9} + {n_{10}D}} \right)}}}{2} \right\rbrack}\end{matrix}{{D = \frac{2G}{{- F} - \sqrt{\left( {F^{2} - {4\; {EG}}} \right)}}},{\beta = \left( {P_{s}/P^{\star}} \right)^{1/4}}}{{E = {\beta^{2} + {n_{3}\beta} + n_{6}}},{F = {{n_{1}\beta^{2}} + {n_{4}\beta} + n_{7}}},{G = {{n_{2}\beta^{2}} + {n_{5}\beta} + n_{8}}}}} & (9)\end{matrix}$

where G is the formation shear modulus in pascals, β is the ratio ofrock matrix to bulk compressibility (c_(rr)/c_(t)), T_(s), is thesaturation temperature in degrees Celsius, and η_(n), are empiricallyderived values referenced by “Release on the IAPWS IndustrialFormulation 1997 for the Thermodynamic Properties of Water and Steam,”The International Association for the Properties of Water and Steam.Erlangen, Germany, September 1997.

The displacement field quantities (u_(r), u_(z)) resulting fromexpansion or contraction of the hydrocarbon bearing reservoir may bedetermined by making the following assumptions: 1) occurrence of linearstress-strain relations, and 2) uniform deformation properties. Based onthese assumptions, the poroelastic stress-strain relation may be givenaccording to Equation 10 as:

$\begin{matrix}{\sigma_{ij} = {{2{G\left\lbrack {e_{ij} + {\frac{v}{1 - {2v}}e\; \delta_{ij}}} \right\rbrack}} - {\left( {1 - \beta} \right)P\; \delta_{ij}}}} & (10)\end{matrix}$

where σ_(ij) is a stress component related to bulk stress system inpascals, e_(ij) maybe a strain component related to bulk stress system,e is Σe_(ij) dilation of the bulk material and δ_(ij) is the Kroneckerdelta.

On the basis of thermoelastic theory, an analogy may be drawn betweenexpansion or contraction due to pressure and temperature having similareffects on the bulk stress-strain system. This analogy may result in athermoelastic stress-strain relationship given by Equation 11:

$\begin{matrix}{\sigma_{ij} = {{2{G\left\lbrack {e_{ij} + {\frac{v}{1 - {2v}}e\; \delta_{ij}}} \right\rbrack}} - {2G\frac{1 + v}{1 - {2v}}\alpha \; T\; \delta_{ij}}}} & (11)\end{matrix}$

For the poroelastic and thermoelastic cases, the transformations may bedenoted by Equations 12 and 13:

$\begin{matrix}{\left. {{c_{b}\left( {1 - \beta} \right)}P_{res}}\Rightarrow{3\alpha \; T_{res}} \right.,{c_{b} = \frac{3\left( {1 - {2v}} \right)}{2{G\left( {1 + v} \right)}}}} & (12) \\{{c_{m,P} = \frac{\left( {1 - \beta} \right)\left( {1 - {2v}} \right)}{2{G\left( {1 - v} \right)}}},{c_{m,T} = {\frac{\left( {1 + v} \right)}{\left( {1 - v} \right)}\alpha}}} & (13)\end{matrix}$

where c_(m) is the uniaxial compaction coefficient.

The stress distribution σ_(ij) should satisfy internal equilibriumconditions, and if the gravity stress field remains almost unaffected,then the equilibrium conditions should follow as σ_(ijj)=0. In the caseof steam or hot water injection, for example, the interest may lie withchanges in strain and stress caused by local increases in excess porefluid pressure above initial reservoir conditions.

The displacement field quantities (u_(r), u_(z)) around a circulardisk-shaped reservoir (e.g., the single well axisymmetric condition) maybe evaluated in terms of Lipschitz-Hankel type integrals I(a,b;d) or themodified type Ĵ_(mnp) involving Bessel functions J_(a,b,m,n) given byEquation 1. If the modified type Ĵ_(mnp) involving Bessel functionsJ_(a,b,m,n) are adopted, then the displacement field quantities may bewritten as follows:

$\begin{matrix}{{{\overset{\sim}{u}}_{r} = \left\lbrack {{{\overset{\_}{J}}_{110}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)} + {\left( {3 - {4v}} \right){{\overset{\_}{J}}_{110}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} - {2z\; {{\overset{\_}{J}}_{111}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}} \right\rbrack}{u_{r} = {{\overset{\sim}{u}}_{r}\left\lbrack {\frac{c_{m_{P,T}}}{2}{h\left( {{{\pm \Delta}\; P},T} \right)}} \right\rbrack}}} & (14) \\{{{\overset{\sim}{u}}_{z} = \left\lbrack {{{sgn}\; \frac{z - H}{R_{P,T}}{{\overset{\_}{J}}_{100}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} - {\left( {3 - {4v}} \right){{\overset{\_}{J}}_{100}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} - {2z\; {{\overset{\_}{J}}_{101}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}} \right\rbrack}{u_{z} = {{\overset{\sim}{u}}_{z}\left\lbrack {\frac{c_{m_{P,T}}}{2}{h\left( {{{\pm \Delta}\; P},T} \right)}} \right\rbrack}}} & (15)\end{matrix}$

where

$\begin{matrix}{{\overset{\_}{J}}_{mnp} = {\left\{ {{sgn}\left( {\zeta - \zeta^{\prime}} \right)} \right\}^{({m + n + p})}J_{mnp}}} & (16) \\{J_{110} = {\frac{2}{\pi \; {k(\rho)}^{1/2}}\left\lbrack {{\left( \frac{2 - k^{2}}{2} \right){K(k)}} - {E(k)}} \right\rbrack}} & (17) \\{{K(k)} = {\int_{0}^{\pi/2}\frac{\theta}{\sqrt{1 - {k^{2}\sin^{2}\theta}}}}} & (18) \\{{E(k)} = {\int_{0}^{\pi/2}{\sqrt{1 - {k^{2}\sin^{2}\theta}}{\theta}}}} & (19) \\{J_{111} = {\frac{\xi \; k}{2\; \pi \; (\rho)^{3/2}}\left\lbrack {{\left( \frac{2 - k^{2}}{2k^{\prime 2}} \right){E(k)}} - {K(k)}} \right\rbrack}} & (20) \\\begin{matrix}{J_{100} = {{- \frac{\xi \; {{kK}(k)}}{2{\pi (\rho)}^{1/2}}} - \frac{k^{2}{\xi^{2}\left( {1 - \rho} \right)}{\Pi \left( {h,k} \right)}}{4{{\pi\rho}\left( {1 + \rho} \right)}} + {1\mspace{14mu} {for}\mspace{14mu} \left( {\rho < 1} \right)}}} \\{= {{- \frac{\xi \; {{kK}(k)}}{2\pi}} + {\frac{1}{2}\mspace{14mu} {for}\mspace{14mu} \left( {\rho = 1} \right)}}} \\{= {{- \frac{\xi \; {{kK}(k)}}{2\pi \; (\rho)^{1/2}}} + {\frac{k^{2}{\xi^{2}\left( {1 - \rho} \right)}{\Pi \left( {h,k} \right)}}{4{{\pi\rho}\left( {1 + \rho} \right)}}\mspace{14mu} {for}\mspace{14mu} \left( {\rho > 1} \right)}}}\end{matrix} & (21) \\{{\Pi \left( {h,k} \right)} = {\int_{0}^{\pi/2}\frac{\theta}{\left( {1 - {h\; \sin^{2}\theta}} \right)\sqrt{\left( {1 - {k^{2}\sin^{2}\theta}} \right)}}}} & (22) \\{J_{101} = {\frac{{k^{3}\left( {1 - \rho^{2} - \xi^{2}} \right)}{E(k)}}{8{\pi (\rho)}^{3/2}k^{\prime 2}} + \frac{{kK}(k)}{2{\pi (\rho)}^{1/2}}}} & (23) \\{{k^{2} = \frac{4\; \rho}{\left( {1 + \rho} \right)^{2} + \xi^{2}}}{h = \frac{4\; \rho}{\left( {1 + \rho} \right)^{2}}}{k^{\prime 2} = {1 - k^{2}}}} & (24)\end{matrix}$

For the sake of convenience, the variables ρ, ζ, and ζ′ adoptedbeginning with Equation 14 have been normalized with respect to thepressure, thermal, or fracture radius

and are given by:

ρ = r , ξ = z , ξ ′ = H ( 25 )

where r is the radial coordinate of interest, z is the depth of interestand H is the source burial depth.

Once the displacement field is known from poroelastic and thermoelasticexpansion or contraction effects, the stress field quantities (σ_(rr),σ_(θθ), σ_(zz), σ_(rz)) may be evaluated from the stress-strainrelations (Equations 10-11). By assuming that the problem is radiallysymmetric, the following relationships outside the reservoir may bederived:

$\begin{matrix}\begin{matrix}{\sigma_{rr} = {{2\; G\frac{\partial u_{r}}{\partial r}} + {\lambda \; e}}} \\{\sigma_{\theta\theta} = {{2G\frac{u_{r}}{r}} + {\lambda \; e}}} \\{\sigma_{zz} = {{2\; G\frac{\partial u_{z}}{\partial z}} + {\lambda \; e}}} \\{\sigma_{rz} = {G\left( {\frac{\partial u_{r}}{\partial z} + \frac{\partial u_{z}}{\partial r}} \right)}} \\{e = {\frac{\partial u_{r}}{\partial r} + \frac{u_{r}}{r} + \frac{\partial u_{z}}{\partial z}}}\end{matrix} & (26)\end{matrix}$

However, the solution may be extended to include the reservoir byadopting a piecewise linear approach to pressure or temperaturegradients. That is, the reservoir thickness may be discretized into anumber of infinitely thin disks and the solution may be integrated withrespect to pressure or temperature within each disk. As with thedisplacement field, the stress field quantities may then also be writtenin terms of the modified type Ĵ_(mnp) Lipschitz-Hankel integralsinvolving Bessel functions J_(a,b,m,n) as follows:

$\begin{matrix}{{\overset{\sim}{\sigma}}_{rr} = {\quad{{\left\lbrack \begin{matrix}{{{\overset{\_}{J}}_{101}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)} + {3{{\overset{\_}{J}}_{101}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} - {2z\; {\overset{\_}{J}}_{102}\left( {\rho,{\xi + \xi^{\prime}}} \right)} -} \\{\frac{1}{r}\left\lbrack {{{\overset{\_}{J}}_{110}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)} + {\left( {3 - {4v}} \right){{\overset{\_}{J}}_{110}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} - {2z\; {{\overset{\_}{J}}_{111}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}} \right\rbrack}\end{matrix} \right\rbrack \sigma_{rr}} = {{\overset{\sim}{\sigma}}_{rr}\left\lbrack {{Gc}_{m_{P,T}}R_{P,T}{h\left( {{{\pm \Delta}\; P},T} \right)}} \right\rbrack}}}} & (27) \\{{{\overset{\sim}{\sigma}}_{\theta\theta} = \left\lbrack {{4v\; {{\overset{\_}{J}}_{101}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} + {\frac{1}{r}\begin{bmatrix}{{{\overset{\_}{J}}_{110}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)} +} \\{{\left( {3 - {4v}} \right){{\overset{\_}{J}}_{110}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} -} \\{2z\; {{\overset{\_}{J}}_{111}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}\end{bmatrix}}} \right\rbrack}{\sigma_{\theta\theta} = {{\overset{\sim}{\sigma}}_{\theta\theta}\left\lbrack {{Gc}_{m_{P,T}}R_{P,T}{h\left( {{{\pm \Delta}\; P},T} \right)}} \right\rbrack}}} & (28) \\{{{\overset{\sim}{\sigma}}_{zz} = \left\lfloor {{- {{\overset{\_}{J}}_{101}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} + {{\overset{\_}{J}}_{101}\left( {\rho,{\xi + \xi^{\prime}}} \right)} + {2\; z\; {{\overset{\_}{J}}_{102}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}} \right\rfloor}{\sigma_{zz} = {{\overset{\sim}{\sigma}}_{zz}\left\lbrack {G\; c_{m_{P,T}}R_{P,T}{h\left( {{{\pm \Delta}\; P},T} \right)}} \right\rbrack}}} & (29) \\{{{\overset{\sim}{\sigma}}_{rz} = \left\lbrack {{{sgn}\frac{z - H}{R_{P,T}}{{\overset{\_}{J}}_{111}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} - {{\overset{\_}{J}}_{111}\left( {\rho,{\xi + \xi^{\prime}}} \right)} + {z\; {{\overset{\_}{J}}_{112}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}} \right\rbrack}{\sigma_{rz} = {{\overset{\sim}{\sigma}}_{rz}\left\lbrack {G\; c_{m_{P,T}}R_{P,T}{h\left( {{{\pm \Delta}\; P},T} \right)}} \right\rbrack}}} & (30)\end{matrix}$

where

$\begin{matrix}{J_{102} = {\frac{k^{3}}{8\pi \; {k^{\prime 2}(\rho)}^{3/2}}\begin{Bmatrix}{{\left\lbrack {\frac{k^{4}\left\lbrack {1 - \left( {\rho^{2} + \xi^{2}} \right)^{2}} \right\rbrack}{4\; \rho^{2}k^{\prime 2}} + 3} \right\rbrack {E(k)}} +} \\{\left\lbrack \frac{k^{2}\left( {\rho^{2} + \xi^{2} - 1} \right)}{4\; \rho} \right\rbrack {K(k)}}\end{Bmatrix}}} & (31) \\{J_{112} = {\frac{k}{2\pi \; (\rho)^{3/2}}\begin{Bmatrix}{{{\frac{k^{2}}{4\; \rho \; k^{\prime 2}}\left\lbrack {\frac{k^{4}\xi^{2}}{k^{\prime 2}} - 1 - \rho^{2}} \right\rbrack}{E(k)}} +} \\{\left\lbrack {1 - \frac{k^{2}{\xi^{2}\left( {2 - k^{2}} \right)}}{8\; \rho \; k^{\prime 2}}} \right\rbrack {K(k)}}\end{Bmatrix}}} & (32)\end{matrix}$

The displacement field quantities (u_(r), u_(z) resulting from adisplacement discontinuity or dislocation disk (e.g., fracturing of thehydrocarbon bearing reservoir) may be determined by making the followingassumptions: 1) the dislocation is disc-shaped and propagates at aconstant distance from the free surface of an elastic half-space, 2) thedislocation is driven by a Newtonian fluid flowing from a point source,3) the elastic half-space is pre-stressed with a uniform lateralcompressive stress (σ₀) and an axial stress (ρgh), and 4) there exists anon-local relationship between net pressure P (r,t) and width W (r,t) ofthe fracture. The non-local relationship stated in assumption 4 may beestablished via the superposition of dislocation disks and may be givenby Equation 33 as:

$\begin{matrix}\begin{matrix}{\left\lfloor {u_{z}^{+} - u_{z}^{-}} \right\rfloor = D_{n}} & {{r < R},{z = H}} \\{\left\lbrack {u_{r}^{+} - u_{r}^{-}} \right\rbrack = {\frac{r}{R}D_{s}}} & {{r < R},{z = H}}\end{matrix} & (33)\end{matrix}$

The previous equation may define the normal D_(n), dislocation and shearD_(s), dislocation components in terms of the fracture surface or facedisplacements. By adopting singular-type solutions for the dislocationdisks, the singular integral equations may be represented as:

$\begin{matrix}{{{\int_{0}^{R}{{G_{nn}\left( {\frac{r}{R},{\frac{s}{R};}} \right)}{D_{n}\left( {s,t} \right)}{s}}} + {\int_{0}^{R}{{G_{ns}\left( {\frac{r}{R},{\frac{s}{R};}} \right)}{D_{s}\left( {s,t} \right)}{s}}}} = {{- \frac{R^{2}}{E^{\prime}}}{p\left( {r,t} \right)}}} & (34) \\{{{\int_{0}^{R}{{G_{sn}\left( {\frac{r}{R},{\frac{s}{R};}} \right)}{D_{n}\left( {s,t} \right)}{s}}} + {\int_{0}^{R}{{G_{ss}\left( {\frac{r}{R},{\frac{s}{R};}} \right)}{D_{s}\left( {s,t} \right)}{s}}}} = 0} & (35)\end{matrix}$

Equation 34 may imply that the normal stress across the plane of thefracture is equal to the negative of net pressure, and Equation 35 mayenforce a condition of shear stress equal to zero on the fracture faces.

With the proper assumptions, as previously described, the displacementfield quantities (u_(r), u_(z)) around a circular disk-shaped reservoir(e.g., the single well axisymmetric condition) due to a prismatic orshear dislocation may be evaluated in terms of Lipschitz-Hankel modifiedtype Ĵ_(mnp) integrals given by Equation 1:

$\begin{matrix}{{\overset{\sim}{u}}_{r} = {\quad{{\begin{bmatrix}{{{- \left( {1 - {2v}} \right)}{{\overset{\_}{J}}_{110}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} + {\left( {\xi - \xi^{\prime}} \right){sgn}\; \frac{z - H}{R_{P}}{{\overset{\_}{J}}_{111}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} +} \\{{\left( {1 - {2v}} \right){{\overset{\_}{J}}_{110}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} - \left\lbrack {\xi - {\left( {3 - {4v}} \right)\xi^{\prime}{{\overset{\_}{J}}_{111}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}} \right\rbrack -} \\{2{\xi\xi}^{\prime}{{\overset{\_}{J}}_{112}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}\end{bmatrix}u_{r}} = \frac{{\overset{\sim}{u}}_{r}W_{F}}{4\left( {1 - v} \right)}}}} & (36) \\{{\overset{\sim}{u}}_{z} = {\quad{{\begin{bmatrix}{{2\left( {1 - v} \right){{\overset{\_}{J}}_{100}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} + {\left( {\xi - \xi^{\prime}} \right){sgn}\; \frac{z - H}{R_{P}}{{\overset{\_}{J}}_{101}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} -} \\{{2\left( {1 - v} \right){{\overset{\_}{J}}_{100}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} - \left\lbrack {\xi + {\left( {3 - {4v}} \right)\xi^{\prime}{{\overset{\_}{J}}_{101}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}} \right\rbrack -} \\{2{\xi\xi}^{\prime}{{\overset{\_}{J}}_{102}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}\end{bmatrix}u_{z}} = \frac{{\overset{\sim}{u}}_{z}W_{F}}{4\left( {1 - v} \right)}}}} & (37)\end{matrix}$

If however, the injection pressure falls below the fracture opening orclosing pressure, as defined by Equation 4, it may then be preferablyassumed that the fracture width goes to zero instantaneously.

Beginning with the singular-type integral equations, as formulated bysuperposing dislocation disk singular objects, the stress fieldquantities (σ_(rr), σ_(zz), σ_(rz)) may then also be evaluated in termsof the modified type Ĵ_(mnp) Lipschitz-Hankel integrals involving Besselfunctions J_(a,b,m,n). Prismatic objects may be considered for someembodiments, although shear objects may also be implemented instraightforward fashion. The radial, normal, and shear stress quantitiesmay be stated as:

$\begin{matrix}{{\overset{\sim}{\sigma}}_{rr} = {\quad{{\begin{bmatrix}\begin{matrix}{{- {{\overset{\_}{J}}_{101}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} + {\left( {\xi - \xi^{\prime}} \right){sgn}\frac{z - H}{R_{P}}{{\overset{\_}{J}}_{102}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} +} \\{{\frac{1}{\rho}\left\lbrack {\left( {1 - {2v}} \right){{\overset{\_}{J}}_{110}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} \right\rbrack} -}\end{matrix} \\{{\frac{1}{\rho}\left\lbrack {\left( {\xi - \xi^{\prime}} \right){sgn}\frac{z - H}{R_{P}}{{\overset{\_}{J}}_{111}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} \right\rbrack} +} \\{{{\overset{\_}{J}}_{101}\left( {\rho,{\xi + \xi^{\prime}}} \right)} - {\left( {\xi - {3\; \xi^{\prime}}} \right){{\overset{\_}{J}}_{102}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} - {2\xi \; \xi^{\prime}{{\overset{\_}{J}}_{103}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} -} \\{\frac{1}{\rho}\left\lbrack {\left( {1 - {2v}} \right){{\overset{\_}{J}}_{110}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} \right\rbrack}\end{bmatrix}\sigma_{rr}} = \frac{{\overset{\sim}{\sigma}}_{rr}W_{f}E}{4{R_{P}\left( {1 - v^{2}} \right)}}}}} & (38) \\{{\overset{\sim}{\sigma}}_{zz} = {\quad{{\begin{bmatrix}{{- {{\overset{\_}{J}}_{101}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} - {\left( {\xi - \xi^{\prime}} \right){sgn}\frac{z - H}{R_{P}}{{\overset{\_}{J}}_{102}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} +} \\{{{\overset{\_}{J}}_{101}\left( {\rho,{\xi + \xi^{\prime}}} \right)} + {\left( {\xi + \xi^{\prime}} \right){{\overset{\_}{J}}_{102}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} + {2{\xi\xi}^{\prime}{{\overset{\_}{J}}_{103}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}}\end{bmatrix}\sigma_{zz}} = \frac{{\overset{\sim}{\sigma}}_{zz}W_{f}E}{4{R_{P}\left( {1 - v^{2}} \right)}}}}} & (39) \\{{{\overset{\sim}{\sigma}}_{rz} = \left\lfloor \begin{matrix}{{{- \left( {\xi - \xi^{\prime}} \right)}{{\overset{\_}{J}}_{112}\left( {\rho,{{\xi - \xi^{\prime}}}} \right)}} +} \\{{\left( {\xi - \xi^{\prime}} \right){{\overset{\_}{J}}_{112}\left( {\rho,{\xi + \xi^{\prime}}} \right)}} + {2{\xi\xi}^{\prime}{{\overset{\_}{J}}_{113}\left( {\rho,{\xi + \xi^{\prime}}} \right)}}}\end{matrix} \right\rfloor}{\sigma_{rz} = \frac{{\overset{\sim}{\sigma}}_{rz}W_{f}E}{4{R_{P}\left( {1 - v^{2}} \right)}}}} & (40)\end{matrix}$

where

$\begin{matrix}{J_{103} = {\frac{k^{3}}{8\pi \; {k^{\prime 2}(\rho)}^{3/2}}\begin{Bmatrix}{{\begin{bmatrix}{{\left( {\frac{k^{4}\left\lbrack {1 - \left( {\rho^{2} + \xi^{2}} \right)^{2}} \right\rbrack}{4\; \rho^{2}k^{\prime 2}} + 3} \right)\left( {\frac{k^{2}{\xi^{2}\left\lbrack {2 - k^{2}} \right\rbrack}}{2\rho \; k^{\prime 2}} - 1} \right)} +} \\{\frac{k^{4}\xi^{2}}{16\rho^{2}k^{\prime \; 2}}\begin{pmatrix}{\frac{2{{k^{2}\left\lbrack {2 - k^{2}} \right\rbrack}\left\lbrack {1 - \left( {\rho^{2} + \xi^{2}} \right)^{2}} \right\rbrack}}{\rho^{2}k^{\prime 2}} +} \\{{17\left\lbrack {\rho^{2} + \xi^{2}} \right\rbrack} - 1}\end{pmatrix}}\end{bmatrix}{E(k)}} +} \\{\left\lbrack {\frac{k^{2}}{4\; \rho}\left\lbrack {1 - \rho^{2} - {2\; {\xi^{2}\left( \frac{k^{4}\left\lbrack {1 - \left( {\rho^{2} + \xi^{2}} \right)^{2}} \right\rbrack}{4\rho^{2}k^{\prime 2}} \right)}} + 3} \right\rbrack} \right\rbrack {K(k)}}\end{Bmatrix}}} & (41) \\{J_{113} = {\frac{k^{3}\xi}{16\; {\pi (\rho)}^{5/2}k^{\prime 4}}\begin{Bmatrix}{{\left\lbrack {\frac{k^{2}{\xi^{2}\left( {2 - k^{2}} \right)}\left( {{8k^{4}} + {3k^{\prime 2}}} \right)}{4\rho \; k^{\prime 2}} - {6\left( {1 - {k^{2}k^{\prime 2}}} \right)}} \right\rbrack {E(k)}} +} \\{\left\lbrack {{\left\lbrack \frac{k^{2}}{2\rho} \right\rbrack 3{k^{\prime 2}\left( {1 + \rho^{2}} \right)}} - {2k^{4}\xi^{2}}} \right\rbrack {K(k)}}\end{Bmatrix}}} & (42)\end{matrix}$

The single-well analysis for the effective displacement, stress,temperature and excess pore pressure field quantities may be extended toa generalized multi-well method through the principles governingsuperposition. The constitutive effects may first be decomposed and thenet displacement and stress field quantities may be calculated. Bypropagating the single-layered half-space solutions via an appropriatepropagation method (e.g., the Thomson-Haskell Propagator Matrix Method)and determining the n-layered solution for a single well, additionalsuperposition of multiple single-well solutions may lead to thegeneralized n-well field-scale solution.

The displacement measurements from tilt arrays or remote sensing may beused to further constrain or improve the layered earth model and layerproperties. Field-wide surveillance methods may include real-timesurveying of earth surface and subsurface displacements via tilt arrays.Another such method may employ remote sensing capabilities (e.g.,Interferometric Synthetic Aperture Radar (InSAR), Light Detection andRanging (LiDAR), Global Positioning System (GPS)) to periodically surveyearth surface displacements. For some embodiments it may be desirable tointegrate field-wide surveillance methods with earth stress analysismethods as part of a calibration scheme and to enable rapid forward orinverse modeling capabilities. The displacement measurements may be usedto further constrain or improve the model and/or properties on the basisof minimizing the square of the error between the measured surfacedisplacements and the displacement field quantity predicted by the earthstress analysis method. Alternatively, a self-calibration or “teach”mode may be introduced into the method whereby the earth model layeringscheme and well log derived layer properties may be iteratively variedbetween practical upper and lower limits until the square of the errorbetween measured surface displacements and the calculated fielddisplacement quantity at z=0 is reduced.

Vertical profiling of μ-seismic events may also be integrated as part ofthe calibration scheme. Vertical profiling of μ-seismic events may beimplemented in either a forward or inverse modeling mode to furtherconstrain calculated displacement and stress field quantities apart fromsurface displacement matching. In the inverse modeling scenario,information from the active or passive monitoring of events (e.g.,source dimension, source magnitude, source location, and elastic strainenergy release) may be used to determine the time-dependent change(e.g., damage or softening) in layer elastic properties. For the forwardmodeling scenario, the stress dependence of layer elastic and inelasticproperties may be prescribed on the basis of experimental formation testdata (i.e., from uniaxial and triaxial geomechanics testing) andinformation about the characteristics of synthetically generatedμ-seismic events may be calculated. In some cases, additionalconstraints on event characterization may be introduced for the forwardmodeling scenario, for example, due to greater uncertainty in predictingthe evolution of μ-seismicity.

As an example of integrating vertical profiling, FIG. 3 depicts a method300 of vertical profiling passive or active μ-seismic events to forecasttime-dependent changes in rock properties with changes in, for instance,the calculated effective temperature, excess pore pressure, displacementand stress field quantities (i.e., an active damage model). The method300 begins at block 302, and various surveillance information may beread at block 304. The information may include electronic signals fromvarious geophone locations 306 and various geophone sonde 308. A sondemay be any subsurface logging tool that carries electrodes, detectors,and the like into a borehole. At block 310 the occurrence of a suitableevent may be determined. If no such event has occurred, data maycontinue to be collected and analyzed until a suitable event does occur.

If such an event has occurred, the collected data may be digitized atblock 312 and incorporated into a velocity model at block 314. At block316, information may be read from a dipole sonde 318, these waveformsmay be digitized, and the digitized waveforms may also be incorporatedinto a velocity model. The velocity model may be used in conjunctionwith a suitable search algorithm to locate hypocenters 322 in athree-dimensional model. Hypocenters may be thought of as the locationwithin the earth where an event occurs. Based on the waveformcharacteristics, source parameters 320, and hypocenter locations, theevents may then be classified 324 as different event types (e.g.,formation heave and shear, casing failures and continuous μ-seismicradiation, which may be triggered to continuously monitor (CMR-T)).

At block 326 the question of whether the event may be classified as aheave may be determined, and if so, the event may be logged with nofollow up at block 328. The various methods used for detecting andmeasuring earth surface displacement previously discussed may be used indetermining and recording a heave event. If the event is classified suchthat a casing failure is indicated at block 330, then a pressure testmay be conducted at block 332 to check for casing integrity. If theevent is classified as continuous μ-seismic radiation at block 334,automatic and continuous monitoring (Autosim) may be initiated at block336. If continuous monitoring positively indicates an event (CMR-E) atblock 338, then Autosim may be implied on subsequent cycles representedat block 342.

The earth stress analysis consists of numerous variables and by applyingμ-seismic data and/or fieldwide surveillance data, the analysis may beconstrained. Constraining the analysis through an integration scheme mayincrease accuracy and responsiveness. One such viable representation ofan integration scheme is shown in FIG. 4, which illustrates a method 400that may enable correlating the predicted injection and productionrelated field quantities (particularly the stress field quantities) atarbitrary depths with vertical profiling of time-lapse μ-seismicity. Themethod 400 begins by reading data regarding fluid properties 406 atblock 402, rock properties 408, injectant pressure 416, injectant flowrate 418, cycle data 410, and stimulation data 412. At block 420 it maybe determined whether there is a well to analyze. If there is no well toanalyze or if all wells have already been analyzed, the method may stop.If, on the other hand, there are more wells to analyze, the method mayproceed to block 422, where it may be determined if there is a timeincrement of data to analyze.

At block 423 it may be determined whether there is flow, and if the flowrate is at or around zero, then injectant pressure may be determined atblock 424. If pressure exists, then fracture extent, fracture width andthermal extent may be calculated at blocks 432 and 434. The pressure andtemperature gradients may then be evaluated at block 436, and theelastic, half-space solution may be determined at block 438.

If it is determined at block 423 that there is a flow rate, then oil andwater flow data 428, 430 may be read at block 426, and productioncalculations may be performed at blocks 440-444, as described above.Whether or not there is flow, continuous monitoring for seismic eventsmay occur at block 446. If an event is detected, it may be digitized atblock 448; analyzed at blocks 450, 452, and 454; and classified at block456.

Depending on the classification of the event (458, 460, 462), automaticand continuous monitoring (Autosim) may be initiated at block 464 andmay be implied on subsequent cycles at block 466 as described herein inreference to FIG. 3. At block 404, depending on the classification ofthe μ-seismic event, changes in the layer mechanical properties may beoptimized to match the change in energy associated with the μ-seismicevents, and the in situ stress/strain state may then be evaluated orupdated accordingly. In other words, the method 400 may compare thepredicted result with the recorded result and iteratively adjust theparameters of the model to produce a more accurate result.

An illustrative example may be of a steam injection process. With therelevant data 406, 408, 416, 418, 410, 412, 428, 430 collected, themethod 400 may determine that fracture will occur at a certain point. Ifan event is detected 468 and that event is determined to have been afracture, the method 400 may then iteratively alter its calculations(i.e., self-calibrate) so that the calculated fracture point matches theactual fracture point.

FIG. 5 depicts a superposed single-well solution 500 for thedisplacement and stress field quantities due to injection-inducedporoelastic, thermoelastic and dislocation effects. The displacementfield quantities u_(z) 516 and u_(r) 518, resulting from the combinationof said effects are depicted in the superposed graph 510. The verticalaxis 512 represents displacement in meters, while the horizontal axis514 represents the radius in meters. For some embodiments, single-welldisplacement quantities may be calculated and then may be incorporatedinto a larger field model.

The superposed single-well solution, as depicted in FIG. 5, permitscalculation of the displacement and stress field quantities at any depth570 within the subterranean overburden, reservoir 566 and underburden568. This is accomplished by assuming the problem of axisymmetry 562about a single well 564 and in-situ stress conditions favoring theinitiation of hydraulically induced horizontal fractures having aplurality of radii 574, comprising a radius of dislocation 572, radiusof poroelastic expansion or contraction 572 and radius of thermoelasticexpansion or contraction 573. Although adopted herein, these assumptionsmay not necessarily be required for general consideration in determiningthe absolute displacement and stress field quantities. For someembodiments, single-well displacement quantities may be calculated andthen may be incorporated into a larger field model.

A field model may consist of data that may be related to a plurality ofsingle wells. Individual well performance and local displacements may beinfluenced by various factors, such as stresses, acting upon theformation due to other wells operating in the same formation. Throughsuperposition, the analysis of individual wells may be combined to moreaccurately model stresses within the formation, the field and theconditions at individual wells. Field models may predict fielddisplacements and, if actual field displacement is measured, then themodel may be checked for accuracy and adjusted so that it betterpredicts the results of the actual event.

Graph 550 graphs superposed stress 552 on a well in GPa versus radius554 in meters, according to the earth stress analysis techniquesdescribed herein. S_(tt) 556 represents tangential stress, S_(rr) 557represents radial stress, S_(rz), 558 represents shear stress and S_(zz)560 represents vertical stress for an example well. For some embodimentsthe calculation of various stresses may allow increased productivity,while potentially avoiding situations in which stress limits may beexceeded. If stress limits are exceeded, damage to valuable equipmentmay occur along with costly delays.

An Exemplary Hydrocarbon Recovery Control Process

It may also be desirable to control hydrocarbon recovery at a “fieldlevel” to improve multi-well interactions while preventing excessivestress or strain-induced well deformations and mechanical failures. Afield-level control process or system may be a variant (eitherlinearized or nonlinear) of the model predictive control (MPC) process,whereby the future behavior of dependent variables (e.g., well operatingconditions) of the dynamic system (well or field-based) may be predictedaccording to past variations or changes in the independent systemvariables (e.g., subterranean layering, layer elastic properties,present well operating conditions, multi-well injection or productionschemes). An advantage of such a process may be that direct or indirectoperating control feedback, on a per-well basis, may be relied on muchless since the dynamic effects of input variations on well mechanicalintegrity will be mostly known a priori.

In FIG. 6, an MPC scheme 600 for forecasting well operating conditionsis depicted. A model predictive controller 610 may be employed toforecast well operating conditions (e.g., rate 620 and/or pressure 622)based on the present risk of compromising well mechanical integrity.Inputs to the model predictive controller 610 may include the layeredearth model 630 of subterranean lithology, properties of the rock 632within the layer (e.g., elastic properties), energetic fluid properties634, and present fluid rate at time t 636 and/or pressure at time t 638.The outputs may be the forecasted fluid rate (at time t+Δt) 620 and/or aforecasted pressure (at time t+Δt) 622.

Using at least some of these inputs, the model predictive controller 610may uniquely calculate earth displacements 640 and stresses 642 alongselected well profiles and may instantaneously evaluate where thecurrent state of stress lies in relation to a well failure envelope,which depicts axial loads, shear loads and a moment. If the currentstress state lies inside the failure envelope 650, a maximum gain inoutput may be predicted iteratively in an effort to minimize errorbetween failure and current stress. Injection and/or production ratesmay be adjusted based on calculations performed regarding the currentstress state's position inside the failure envelope. These adjustmentshave the potential to increase productivity while reducing the chance ofcostly failures.

If the current stress state falls outside the envelope, the output maybe triggered to a “wait” or “off” state where operations are temporarilyhalted. A wait state may be maintained until the current stress statesreturns to a safe position inside the failure envelope. However, anotherscenario may be when the controller predicts the intersection of thewell failure envelope and generates an alternate scenario stress-strainprediction (ASP) on-line to avoid intersecting the well failure envelopeand triggering a “wait” state. An alert may also be given; allowing anopportunity for an operator to adjust production parameters manually.The ASP should include appropriate constraints on inputs (e.g., boundsfor well operating conditions including number of active wells,subterranean layer elastic properties, and number of layers in the earthmodel representation).

In FIG. 7, a model predictive control scheme 700 that may utilize earthdisplacement measurements to self-calibrate is depicted. A modelpredictive controller 710 may utilize tilt array or remote sensingmeasurements (e.g., Interferometric Synthetic Aperture Radar (InSAR),Light Detection and Ranging (LiDAR), Global Positioning System (GPS)) ofearth surface displacements 720 in real-time to further constrain themodel and self-calibrate (“teach”) prior to forecasting fluid rate (attime t+Δt) 620 and/or forecasting pressure (at time t+Δt) 622. As anexample, if the model-based calculations for earth displacements 640 atground surface match the measured surface earth displacements, the MPCmay continue to forecast. Otherwise, rock property inputs 632 and theearth model layer scheme 630 may be varied autonomously within thespread of reported and/or experimental values (“input constraint”) untila reduced error is achieved between the calculated and measureddisplacements. Once a desired minimum error is reached, the modelpredictive controller 710 may proceed with evaluating the output on thebasis of comparing ASPs for the possible permutations of constrainedlayer elastic properties (e.g., log-based and/or geomechanicstest-based) and layering schemes (e.g., log-based and/or core-based)given current well operating conditions.

In FIG. 8, a model predictive control scheme 800 that may utilize earthdisplacement measurements and vertical profiling of microseismic(μ-seismic) events is depicted. The model predictive controller 810 mayutilize measurements of both earth displacements 820 and verticalμ-seismic profiling 830 simultaneously prior to forecasting welloperating conditions (at time t+Δt). Within the MPC scheme 800, verticalμ-seismic event profiling 830 may be implemented in either a forward orinverse mode in an effort to further constrain forecasting displacement(or strain) and stress calculations aside from surface displacementmatching. In an inverse modeling mode, information from the active orpassive monitoring of events (e.g., source dimension, source magnitude,source location, elastic strain energy release) may be used to predictthe time-dependent change (e.g., damage or softening) in layer elasticproperties (rock property inputs 632).

In a forward modeling mode, the stress dependence of layer elastic andinelastic properties may be prescribed on the basis of experimentalformation test data (e.g., from uniaxial and triaxial geomechanicstesting), and information about the characteristics of syntheticallygenerated μ-seismic events may be predicted. Additional constraints onevent characterization may be required for the forward modeling scenariobecause of greater uncertainty in predicting the evolution ofμ-seismicity.

Well mechanical integrity may be managed by a well-located MPG system.One element of a well-located MPG system may be a physics-based control“engine” for transient analysis of formation effective displacement,stress and excess pore pressure field quantities. A strain or stressbased systematic method for analysis of the multi-well problem throughdecomposition of phenomena governing single-well mechanical response, asdescribed above, may be used in an MPC system.

An Exemplary Computer-Based Model Predictive Control Scheme

In FIG. 9, a computer-based implementation 900, using an interactivegraphical computer software code may be designed to be flexible andinterpretive in its use. The code may consist of a graphical userinterface (GUI) comprising a variety of modules. As an example, themodules may be used interactively to manage data inputs at block 910,construct the layered earth model and layer elastic properties at block920, calculate and view field quantities at block 930, evaluate wellmechanical integrity at block 940, collect field data at block 950 andself-calibrate at block 960. However, other computer-basedinstantiations of the MPC scheme may be implemented. Accordingly, thefollowing description of the example GUI is intended not to beconsidered as limiting, but to be illustrative as an exampleimplementation of the predictive control system or process.

In FIG. 10, an exemplary main GUI 1000 may enable functionality fordefining a project, the type of analysis (i.e., a single-layer ormulti-layer analysis) and the paths for solution data retrieval andstorage. The main GUI may also enable various exemplary modules asfollows: a data management module 1010, a well log calculation module1020, a field quantity calculation module 1030, a field visualizationmodule 1040, a well integrity module 1050 and an MPC self-calibrationmodule 1060. Below is a brief description of the basic functionality ofthese modules.

FIG. 11 illustrates an exemplary data management module 1100 that maydefine the type of operational process to evaluate and the associatedinput formation and injectant properties. An example may be asteam-based thermal recovery process where steam may be injected in amulti-well row-by-row scenario (i.e., a pad configuration), hence steamproperties and multi-pad selection may be permitted. FIG. 12 depicts adata management module 1200 according to another embodiment of theinvention that may permit the user to graphically select the time frameand associated well operating parameters over which a calculation offield quantities is to occur. In a control mode, the time frame may besynchronized with the frequency at which the well operating parametersmay be sampled.

As depicted in FIG. 13, an exemplary well log calculation module 1300,may be enabled when the analysis type is a multi-layered analysis. Thismodule may permit the user to load log files from representative wellsfor use in evaluating the layered earth model and layer elasticproperties. In this example, the stratified earth model is defined fromthe gamma ray log, although multiple logs or composites from multiplelogs (e.g., bulk density, resistivity and sonic) may be used in definingthe earth model. Once the earth model has been stratified via aconvolution of the logs, then layer properties may be calculated usinganalytical relationships or empirical correlations known to thoseskilled in the art.

Once the earth model and associated layer properties have beendetermined, a transient analysis of field quantities based onconstitutive effects and input data may be calculated. FIG. 14 is anexemplary depiction of the field quantity calculation module 1400, whichmay permit interrogation of the single-well solution (eitherporoelastic, thermoelastic, dislocation or some combination thereof) forany particular well. This module may display the transient results forthe fracture and thermal fronts, and displacement and stress componentsat any arbitrary subterranean depth. FIG. 15 depicts another exemplaryform of the field quantity calculation module 1500, which may allow theuser to graphically define the areal extent over which the single-wellsolutions may be superposed to solve for the field quantities. This formmay be flexible enough to compensate for various projection standards(e.g., State Plane Coordinate System (SPCS), 3-degree TransverseMercator (3TM), the Universal Transverse Mercator (UTM), and theGeodetic Reference System (GRS90)) and datum references (e.g., NorthAmerican Datum 1927 (NAD 27), North American Datum 83 (NAD 83), andWorld Geodetic System 1984 (WGS 84)), and this form may also permitcoordinate transformations.

FIG. 16 depicts a field visualization module 1600 that may enable theuser to visualize the spatial variation of field quantities (e.g.,excess pore pressure, displacement, and strain or stress via a “Setup”menu) as a function of subterranean depth. This module may also enablethe user with a capability to simultaneously rotate or scale the 3-Dfield display in any orientation, and it may also be capable ofinterrogating (querying) the superposed solution at some selected pointor coordinate in space for the interpolated result. Since the full-fieldsolution should now be known, the user may be able to evaluate wellmechanical integrity on a “per well” basis.

FIG. 17 depicts a well integrity module 1700 that may be employed forquerying particular components of displacement, strain or stress alongindividual well profiles, which may comprise the geometric description(geometry data file) in conjunction with well completion details (e.g.,completion type, cement type, casing grade and weight). Like the fieldvisualization module 1600, the well integrity module 1700 may havesimilar 3-D display capabilities, but permit the user to interpretdepth-dependent well deformation and strain components (e.g., flexuralstrains) either directly inferred or calculated.

The well integrity module 1700 may also have on-line warning or alarmfunctionality whereby the user is notified in the event that a single ormultiple wells have met certain mechanical integrity criteria (e.g., thebuckling limit or the shear-slip limit). FIG. 18 depicts an adoption1800 of a shear-slip limit approach as a constraint on integrity for theexample depicted in FIG. 17. That is, if a certain magnitude and“interaction” of the radial and vertical component of displacement isdetected along any well profile, then a warning and/or alarm stateshould be displayed to the user. In control mode, this state could leadto a change in the control state (e.g., a gain reduction) for theinjection or production rates. In an MPC-based configuration, theaforementioned concept may easily be envisioned to adopt quantitativerisk analysis measures for constraints on prediction of alternateinjection or production operating scenarios.

FIG. 19 depicts an implementation of a model predictive control (MPC)self-calibration module 1900, which may enable the functionality forprocessing remote sensing measurements (e.g., Interferometric SyntheticAperture Radar (InSAR), Light Detection and Ranging (LiDAR), and GlobalPositioning System (GPS)). As an example, the module may provide thecapability to create interferograms from 2-pass, single-look complex(SLC) synthetic aperture radar (SAR) images. From the interferograms ofthe two SLC SAR images, time-dependent vertical surface displacements ofthe ground surface may be employed in an MPC-based control configurationfor predictive model self-calibration. Measured surface displacementsand, consequently, surface strains or stresses may also be compared withthe stress inversion of microseismic events for auto-correction ofpredicted spatial distribution of stress.

An Exemplary Fluid Injection Process

A variant (either linearized or nonlinear) of the model predictivecontrol (MPC) process may also be used for controlling fluid injectionparameters in an effort to improve well interactions and controlhydrofracture geometries. For example, the future behavior of dependentvariables (e.g., well operating conditions) of the well or field-baseddynamic system may be predicted according to the past variations orchanges in the independent system variables (e.g., subterraneanlayering, layer elastic properties, present well operating conditions,multi-well injection or production schemes). An advantage of thisprocess may be that direct or indirect operating control feedback, on aper-well basis, may be relied on much less since at least some of thedynamic effects of input variations on well mechanical integrity maylikely be known a priori according to the earth stress analysis methoddescribed herein. While those skilled in the art will recognize thatvarious fluids may be used in injection operations (e.g., steam, carbondioxide, acid, and natural gas), water will be used henceforth as anexemplary injectant.

FIG. 20 depicts an MPC scheme 2000 to forecast well operating conditionsrelative to injection constraints. A model predictive controller 2010may be employed to forecast future injector operating conditions 2030,which may include pressure and/or rate relative to current injectionperformance. The model predictive controller inputs may include thelayered earth model 630 of subterranean lithology, rock properties 632(e.g., layer elastic properties, layer strength properties), energeticfluid properties 634, and present fluid rate 636 at time t and/orpressure 638 at time t. The outputs may be the forecasted fluid rate attime t+Δt 620 and/or pressure at time t+Δt 622.

Using the inputs, the model predictive controller may uniquely calculatethe convolution of fracture growth and adapt tilt array or remotesensing (e.g., InSAR, LiDAR and GPS) of earth displacement measurements2020 in real-time to constrain the calculation. Calculated fracturegrowth may then be compared to a target fracture extent, and injectionparameters, such as gain, may be adjusted to reduce the error.

FIG. 21 depicts an MPC model predictive control scheme 2100, which mayemploy earth displacement measurements 2120 to self-calibrate. A modelpredictive controller 2110 may calculate earth displacements 640 andstresses 642 and may predict a likely injection mode 2130 (e.g., matrixinjection, fracturing or fluidization). If the mode is matrix injection,a maximum gain may be specified until a change in mode is detected. Ifthe mode is fracturing, the output gain may be predicted iterativelyaccording to a maximum constraint on the calculated convolution offracture growth. If the mode is fluidization, the model predictivecontroller 2110 may utilize the measurements of earth surfacedisplacements 2120 to predict the radial and vertical extent of thenear-well disturbance. If the extent of the disturbance lies insidebounding strata, the model predictive controller 2110 may continue toforecast a gain in output. However, if the extent approaches boundingstrata and pressure predicted within the extent exceeds the strength ofthe strata, the output may be triggered to a “wait” or “off” state.

In FIG. 22, an MPC scheme 2200 is depicted which may employ earthdisplacement measurements 2210 and vertical profiling of microseismic(μ-seismic) events 2240. The model predictive controller 2210 mayutilize measurements of earth displacements 2210 and vertical profilingof injection-induced μ-seismic events 2240 simultaneously prior toforecasting at time t+Δt. Within the model predictive controller 2210,vertical μ-seismic event profiling may be implemented in an effort tocalibrate the forecasted mode of injection and may further constrain thecalculated convolution of fracture growth (if in fracturing mode) andextent of any near-well disturbance (if in fluidization mode).Information from the monitoring of events (e.g., source dimension,magnitude, location and elastic strain energy change) may be used todetermine at t+Δt when a change in injection mode is expected andwhether or not the t+Δt mode turns out to be aseismic (i.e., the rock ismore fluidic rather than intact).

Water injection may be managed by a well-located model predictivecontrol (MPG) system. One element of a well-located MPG system is aphysics-based control “engine” for transient analysis of formationeffective displacement, stress and excess pore pressure fieldquantities. A strain or stress based systematic method for analysis ofthe multi-well problem through decomposition of phenomena governingsingle-well mechanical response, as previously described, may be used inan MPC system to control water injection.

Those skilled in the art should understand that the preferred embodimentherein discloses a control system or process that is preferablyimplemented for field-wide management of water injection using asuitably programmed digital computer. Such persons could develop acomputer software and hardware implementation of the invention based onthe methods described herein for management and control of earth stress.

While the foregoing is directed to embodiments of the present invention,other and further embodiments of the invention may be devised withoutdeparting from the basic scope thereof, and the scope thereof isdetermined by the claims that follow.

1. A method for managing an impact, on an earth formation, ofhydrocarbon recovery operations from at least one well formed in theearth formation and preventing harm to the well, comprising: a)generating at least first and second sets of equations to modelcontributions to the impact of the operations due to at least first andsecond physical processes associated with the operations, whereingenerating the at least first and second sets of equations comprisesanalyzing at least one of historical data related to the hydrocarbonrecovery operations, energetic fluid properties, rates of fluids beingproduced or injected during the operations, pressures of fluids beingproduced or injected during the operations, layer elastic properties,and layered earth models of subterranean lithology; b) obtainingsolutions to the first and second sets of equations to determinecontributions to the impact of the operations due to the first andsecond physical processes; c) combining the solutions to the first andsecond sets of equations to determine the impact of the operations onthe earth formation; and d) adjusting the hydrocarbon recoveryoperations of the well based on the combined solutions.
 2. The method ofclaim 1, further comprising: dividing the well into a plurality oflayers; conducting steps a-c for each of the plurality of layers togenerate a plurality of combined solutions for the layers; superposingthe plurality of combined solutions to determine the impact of theoperations at the well on the earth formation; and adjusting thehydrocarbon recovery operations of the well based on the superposedsolutions.
 3. The method of claim 1, further comprising: repeating stepsa-c for a plurality of wells to generate a plurality of combinedsolutions for the wells; superposing the plurality of combined solutionsto determine a field-level impact of the operations on the earthformation; and adjusting the hydrocarbon recovery operations of theplurality of wells based on the superposed solutions.
 4. The method ofclaim 1, further comprising: e) dividing the well into a plurality oflayers; f) conducting steps a-c for each of the plurality of layers togenerate a plurality of combined solutions for the layers; g)superposing the plurality of combined solutions for the layers todetermine the impact of the operations at the well on the earthformation; h) repeating steps e-g for a plurality of wells to generate aplurality of combined solutions for the wells; i) superposing theplurality of combined solutions for the wells to determine a field-levelimpact of the operations on the earth formation; and j) adjusting thehydrocarbon recovery operations of the plurality of wells based on thesuperposed solutions for the wells.
 5. The method of claim 1, whereinthe first and second sets of equations are generated based, at least inpart, on historical data related to the hydrocarbon recovery operations.6. The method of claim 5, wherein the historical data comprises at leastone of variations in subterranean layering, variations in layer elasticproperties, variations in present well operating conditions, variationsin multi-well injection schemes, and variations in production schemes.7. The method of claim 1, wherein generating the first and second setsof equations comprises analyzing at least one of a layered earth modelof subterranean lithology, layer elastic properties, energetic fluidproperties, a rate of a fluid being produced or injected during theoperations, and a pressure of the fluid.
 8. The method of claim 1,further comprising determining a well failure envelope based on acritical set of stresses under which the well will fail.
 9. The methodof claim 8, wherein the critical set of stresses comprises one or moreof an axial loading, a shear loading, and a moment.
 10. The method ofclaim 8, further comprising determining where a current stress statelies in relation to the well failure envelope.
 11. The method of claim10, further comprising iteratively predicting a maximum gain andreducing an error between a failure state and the current stress stateif the current stress state lies within the well failure envelope. 12.The method of claim 10, further comprising temporarily halting thehydrocarbon recovery operations if the current stress state is outsideor near a boundary of the well failure envelope.
 13. The method of claim10, further comprising predicting an intersection of a projected stressstate and a boundary of the well failure envelope and generating analternate scenario to avoid the intersection.
 14. The method of claim 1,further comprising using earth surface displacement measurements toconstrain a well model based on the combined solutions to the first andsecond sets of equations.
 15. The method of claim 14, wherein usingearth surface displacement measurements to constrain and calibrate thewell model comprises: iteratively comparing a predicted earthdisplacement to an actual earth displacement; adjusting rock propertydata to update the model; adjusting earth model layer scheme data toupdate the model; and halting the iterative comparisons when thepredicted earth displacement is sufficiently similar to the actual earthdisplacement.
 16. The method of claim 14, wherein the earth surfacedisplacement measurements are from one or more remote sensing devicescomprising at least one of Interferometric Synthetic Aperture Radar(InSAR), Light Detection and Ranging (LiDAR), and Global PositioningSystem (GPS) devices.
 17. The method of claim 1, further comprisingforecasting at least one of a future fluid rate and a future fluidpressure of a fluid being produced or injected during the operations.18. A computer-readable storage medium containing a program for managingan impact, on an earth formation, of hydrocarbon recovery operationsfrom at least one well formed in the earth formation and preventing harmto the well, which when executed performs operations comprising:generating at least first and second sets of equations to modelcontributions to the impact of the operations due to at least first andsecond physical processes associated with the operations, whereingenerating the at least first and second sets of equations comprisesanalyzing at least one of historical data related to the hydrocarbonrecovery operations, energetic fluid properties, rates of fluids beingproduced or injected during the operations, pressures of fluids beingproduced or injected during the operations, layer elastic properties,and layered earth models of subterranean lithology; obtaining solutionsto the first and second sets of equations to determine contributions tothe impact of the operations due to the first and second physicalprocesses; combining the solutions to the first and second sets ofequations to determine the impact of the operations on the earthformation; and adjusting the hydrocarbon recovery operations at the wellbased on the combined solutions.
 19. A system for managing an impact, onan earth formation, of hydrocarbon recovery operations from at least onewell formed in the earth formation and preventing harm to the well,comprising: a processing unit configured to (a) generate at least firstand second sets of equations to model contributions to the impact of theoperations due to at least first and second physical processesassociated with the operations, wherein generating the at least firstand second sets of equations comprises analyzing at least one ofhistorical data related to the hydrocarbon recovery operations,energetic fluid properties, rates of fluids being produced or injectedduring the operations, pressures of fluids being produced or injectedduring the operations, layer elastic properties, and layered earthmodels of subterranean lithology, (b) obtain solutions to the first andsecond sets of equations to determine contributions to the impact of theoperations due to the first and second physical processes, (c) combinethe solutions to the first and second sets of equations to determine theimpact of the operations on the earth formation, and (d) adjust thehydrocarbon recovery operations at the well based on the combinedsolutions.
 20. The system of claim 19, wherein the processing unit isfurther configured to repeat steps a-c for a plurality of wells togenerate a plurality of combined solutions for the wells, superpose theplurality of combined solutions to determine a field-level impact of theoperations on the earth formation, and adjust the hydrocarbon recoveryoperations of the plurality of wells based on the superposed solutions.